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Abstract 

Coherent angular rotation of epithelial cells is thought to contribute to many vital physiological processes including 
tissue morphogenesis and glandular formation. However, factors regulating this motion, and the implications of 
this motion if perturbed remain incompletely understood. In the current study, we address these questions using 
a cell-center based model in which cells are polarized, motile, and interact with the neighboring cells via harmonic 
forces. We demonstrate that, a simple evolution rule in which the polarization of any cell tends to orient with its 
velocity vector can induce coherent motion in geometrically confined environments. In addition to recapitulating 
coherent rotational motion observed in experiments, our results also show the presence of radial movements and tissue 
behavior that can vary between solid-like and fluid-like. We show that the pattern of coherent motion is dictated 
by the combination of different physical parameters including number density, cell motility, system size, bulk cell 
stiffness and stiffness of cell-cell adhesions. We further observe that, perturbations in the form of cell division can 
induce a reversal in the direction of motion when cell division occurs synchronously. Moreover, when the confinement 
is removed, we see that the existing coherent motion leads to cell scattering, with bulk cell stiffness and stiffness 
of cell-cell contacts dictating the invasion pattern. In summary, our study provides an in-depth understanding of 
the origin of coherent rotation in confined tissues, and extracts useful insights into the influence of various physical 
parameters on the pattern of such movements. 


Introduction 

Collective cell migration is central to both physiological processes such as morphogenesis and wound healing, and 
pathological processes like cancer invasion w- Epithelial and endothelial cells collectively migrate in intricate pat¬ 
terns within a tissue by virtue of their adhesion to their neighboring cells and to the extracellular matrix (ECM) [7, 8]. 
Further, on 2D confined geometries, these cells exhibit coherent angular movement (CAM) [9UI2]. Interestingly, such 
coordinated movements have also been documented in various in vivo processes including egg chamber elongation 
in Drosophila , ommatidial rotation in Drosophila and development of spherical mammary acini [T3HT8] . In addition 
to these types of tissues, such large scale rotations are also observed in confined dictyostelium colonies and bacterial 
suspensions DSEDl- Moreover, non-living, active materials such as vibrated, granular materials also exhibit sponta¬ 
neous CAM when confined ED- Thus, large scale rotational movements under confinement are ubiquitous in ‘active 
systems’ - both non-living and living. 
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Active systems have been modeled using a variety of approaches ranging from discrete, self-propelled particle 
modeling (SPM) to active hydrodynamical theories [22j[23]. Of special interest are theories, which involve discrete 
or continuum elements with self-propulsion, and are successfully used to describe collective motion in epithelia [HI, 
[26h29l[83l[93] . The common thread connecting these diverse modeling attempts is the presence, in some form, of self 
propulsion velocity Vo originating from actin polymerization and polarization p for the active elements, in addition 
to the elastic and viscous interactions of the elements with their surrounding constituents. The polarization p is a 
coarse-grained representation of front-rear asymmetry of a migrating cell resulting from various factors, e.g., Rho 
GTPase gradient [30] and position of centrosome in relation to the nucleus [HUES- A SPM-based cellular Potts 
model has successfully replicated the existence of CAM in confined epithelia [9]. Similarly, a recent study has also 
demonstrated that a particle based model for confined epithelia, where cells are represented as self-propelled points 
connected to their neighbors with elastic springs, also gives rise to CAM [93] • A distinct, but related, formalism that 
utilizes dissipative particle dynamics for sub-cellular components, has been used to demonstrate spontaneous rotation 
of two tightly connected, and confined cells [33]. Recently, Camley et al. used a phase-field method for studying 
the emergence of coherent rotation in a pair of cells confined to adhesive micropatterns, and demonstrated that even 
subtle differences in cell polarity mechanisms (via contact inhibition, neighbor alignment, velocity alignment, etc) 
greatly influence the pattern of collective movement m ■ 

Though SPM has also been utilized to represent stable vortex formation in confined bacterial suspensions and 
driven granular media [20]|2T], the mechanism of vortex formation relies on hydrodynamic coupling between the 
active particles through the surrounding media. This is distinct from the collective behavior in epithelia, in which 
the tight inter-cellular contacts play a crucial role in tissue movements m\ ■ 

Despite the presence of several models addressing CAM in epithelia, a number of crucial questions still remain 
unanswered. For example, there is no simple understanding as to why CAM spontaneously emerges in such systems. 
Additionally, the diversity of experimental findings from similar experiments with near-identical setup [9] EH [83] 
raises the possibility that small perturbations in physical parameters associated with the experiments are likely to 
influence the various hydrodynamic modes exhibited by cells, and hence motivates the necessity of identifying some 
of the critical physical parameters. Though the confinement provided for in vitro cultures arises quite naturally 
on the micropatterned geometries, the confinement for epithelia in vivo comes from being embedded in a larger 
tissue [351136] . What role the nature of confinement plays on the emergence and sustenance of CAM is another 
issue that needs addressing. Similarly, the influence of confinement geometry on CAM is also not clear, and is 
particularly relevant to several in vivo situations. For example, the annular geometry is the simplest non-convex 
geometry that is relevant in understanding CAM in biological lumens m ■ Though CAM on this geometry has been 
addressed before m , the role of cell cohesivity on the nature of CAM still remains unaddressed. Finally, the roles of 
internal and external perturbations mimicking various in vivo processes, in the form of cell division and the loss of 
confinement on CAM also remain unknown. In this paper, we employ SPM using cell center representation of cells 
to computationally answer these questions [261193] . Additionally, using simple calculations we also provide analytical 
insights to get a better understanding of CAM in epithelia. We show that the nature of coherent motion can be 
both solid-like or fluid-like, and is dictated by the combination of different physical parameters including number 
density, cell motility, system size, bulk cell stiffness and stiffness of cell-cell adhesions. Our results also predict that 
synchronous cell division can lead to change in the direction of rotation for CAM in annular geometries. Finally, we 
show that depending on the properties of cells and cell-cell adhesions, CAM leads to different patterns of invasion 
when the confinement is removed. Collectively, our results illustrate the influence of cell density, system size, cell 
motility, cell division, and stiffness of cell and cell-cell adhesions in regulating CAM. 


Computational Model 

An epithelial sheet is comprised of a group of cells that are connected to each other via cadherin bonds to form a 
monolayer. Many experimental observations have demonstrated that cells in this network are persistently motile, 
and upon reaching a critical density show collective migration behavior [38]. Presence of front-rear polarity axis is 
known to be essential for migrating cells. This polarity axis manifests in migrating cells in different forms like: (i) 
increased actin activity in the front and formation of actin structures such as lamellipodia, (ii) localization of the 
microtubule organizing center (MTOC) at the front of the nucleus with microtubule growth towards the leading 
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edge, (iii) gradients in cell-ECM adhesion, and (iv) establishment of front-rear gradients in the activity of GTPases 
such as Rac/Cdc42 [39]. Cell polarity is actively maintained and constantly steered by complex mechano-chemical 
processes governed by cell-cell and cell-ECM interactions [40][4I|. A surprisingly simple upshot of these complex 
processes in terms of mechanical observables is that, in epithelial sheets such as MDCK tissue, the polarization of 
constituent cells is closely oriented with the principal direction of stress as well as with their average velocity @203]. 
Keeping these experimental observations in mind, we have utilized a simple model to explore how mechano-chemical 
properties of individual cells impact their collective behavior in confined epithelial sheets. 

For modeling the collective mechanics of cells, we have adopted a ‘cell center-based mechanics model’ with cells 
represented as discrete points at their center of mass 26. 1451194] . As shown in Fig. [1] the whole epithelial tissue is 
represented as a continuous sheet with cell-cell cadherin junctions represented by simple harmonic springs [261146] , 
Each cell is assumed to exert an attractive or repulsive force on its neighboring cells depending on the relative 
deformation of springs with respect to their undeformed length, ao and stiffness, k. The force acting on any cell at 
any time, £, is the sum of the contributions of all the connecting neighbors. Thus, if ri represents the position of i th 
cell, the net force exerted on that cell by neighbors (m, say) is given by 

Fi = ^2 k (\ r i ~ r *l “ a o) e ij (!) 

j'G neighbor 


where, represents the unit vector along the direction connecting the i th cell with its j th neighbor. 

Depending on the relative deformation of springs with respect to the natural length, the interaction potential can 
either be tensile or compressive. In order to avoid force transfer between distant neighbors, it is assumed that when 
the deformation of spring is greater than a threshold, d max , no force transfer occurs between those two cells. For all 
our simulations, we took the value of d max equal to 1.3 ao [26] . Thus the value of spring stiffness for the entire range 
of deformation can be written as: 


k = < 


0 , 

kt: 

kn i 


If (IT? ^i I do) ^ ^max- 

if 0 (| v*j v 1 1 ao) ^ d max . 

if (|rj - ri \ - a 0 ) < 0. 


( 2 ) 


In the above expression, k c and k t represent the bulk cell stiffness and the stiffness of cell-cell adhesions (or cohesivity), 
respectively. Figure QJe) illustrates the attractive/repulsive force acting on each cell. The cells are allowed to 
exchange their neighbors, which are obtained by repeated Delaunay triangulation |93] l94]. For a given set of cell 
centers, Delaunay triangulation provides a connectivity for cells that produces the least number of distorted triangles, 
i.e., triangles with least shear strain. Delaunay triangulations are dual to Voronoi tessellations (Fig. [T| (b), (c)) and 
the Voronoi polygon for a given cell center can be modeled to be the cell itself (see Materials and Methods). 
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Figure 1. A schematic of cell center model depicting the arrangement of cells and the forces acting 
on them, (a) A 2-D monolayer of epithelial cells, confined inside a circular geometry is considered with cells 
represented as points at their center, (b) Delaunay triangulation (blue) has been used to model cell - cell 
connectivity, which finds the nearest neighbors of each point and form the connectivity array accordingly. Because 
of the greater clarity it affords and better connection with the experimental geometry, Voronoi tessellation 
(topological dual of Delaunay triangulation) is used for visualization of cells, (c) When two originally connected 
cells move apart and form new neighbors, the connectivity of the system is updated using Delaunay triangulation. 
This connectivity update automatically takes T1 transitions into account, (d) Enlarged view of a representative cell 
i, along with its connection to neighboring cells. The position vector of this cell center is denoted by r* and position 
vector of its j th neighbor is denoted by rj. The blue arrow indicates the force, F^ acting between cells i and j. The 
total force acting on i th cell is the sum of the contributions from all the connecting neighbors, (e) The interaction 
between two adjacent cells is either compressive or tensile, depending upon the relative deformation of connecting 
spring with respect to its undeformed length, ao- Here compressive and tensile stiffness of each spring is 
represented by k c and kt , respectively. While k c mimics the bulk cell stiffness, kt mimics cell-cell cohesivity. It is 
assumed that if the deformation of any spring is greater than d max , the cell-cell connection is broken and there is no 
force transfer between these two cells, (f) Force acting on each cell is resolved along anti-parallel (F\\) and 
perpendicular (F±) to the direction of the cell’s polarization(p). Here v denotes the velocity vector on each particle, 
(g) Velocity profile in the direction of polarization as a function of F\\. 


In our model, cells are assumed to act as self propelled active particles [22], with their inherent motility (uo) 
representing the speed with which they move in the absence of any external force. The preferential direction of cell’s 
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motion (i.e., polarization) is represented by the vector p*, which is a coarse-grained representation of the front-rear 
polarization in a motile cell [39]. As cells move over a viscous substrate with mobility /i, the drag force acting in the 
opposite direction of motion balances the internal forces. If ri is the position vector of i th cell, its velocity at time t 
can be written as: 

dv' 

Vi = —j£ = Vo Pi + fiFi (3) 

Similar to the procedure followed elsewhere m and as motivated earlier, we assume that the cell’s polarization 
vector tends to orient with its velocity vector as per the following equation: 


dpi 

dt 


£(Pi x Vi.e z )pj- 


(4) 


where Vi is the unit velocity vector and e z is the unit vector perpendicular to the plane. The parameter £ represents 
the polarization coordination constant determining the tendency of cell’s polarization to rotate and align with the 
velocity vector. We do not account for noise in our simulations m as we are primarily interested in a mean-field 
understanding of CAM, and noise is known to typically increase fluctuations in the system [931 . 


Numerical estimates of parameters used in the study 

Before employing our model for any qualitative predictions, it is essential to estimate the real values of various 
parameters used in the study. Consistent with previous studies mmm, 20 /i m and 20 /rm/hr were taken as cell 
length and cell speed, respectively. Unless specified, for all the simulations, the non-dimensional value of ao and 
vo were taken as 1. Assuming a substrate drag coefficient (£) of 100 pN hr//im 3 [27], the value of mobility /i was 
calculated to be l/(£ag) = 2.5 x 10 -5 /inl/pN hr. Length, time and force are expressed in units of ao = 20 /im, 
To = clo/vq = 1 hr and /o = vq/ji = 8 x 10 5 pN, respectively. The non-dimensional value of mobility (fi) was taken 
as 1 for all the simulations unless otherwise specified. The value of stiffness of cell-cell connection is given by the 
expression k = E h/(2y/3 (1 — v)) [79] where h, E and v represent the height, Young’s modulus and Poisson’s ratio 
of the cell, respectively (see supplementary information (SI) Text SI for derivation). Assuming values of E = 10 kPa, 
h = 5 fi m and v = 0.5 [43], the value of stiffness was estimated to be k « 0.03 N/m. In our simulations, we have 
used the non-dimensional value of stiffness (k) in the range of 1 — 10. For this range, the real value of stiffness was 
calculated as k m kvo/fiao yielding a value of k = 0.04 — 0.4 N/m, which is close to the actual value of cell-cell 
stiffness derived above. Due to uncertainty in the value of £, the non-dimensional value of k can indeed have a some 
variability. The outer radius of substrate was taken as 100 fi m for all the simulations [9], unless stated otherwise. 
For annular geometry, the inner radius was taken as 70 fi m. The number of cells were varied between 100 — 170 for 
various simulations, yielding an average cell density in the range of 3000 - 6000 cells/mm 2 which closely matches 
with previous experimental studies [9]. 


Results 

Coherent rotation of cells confined in circular geometry 

Various theoretical studies modeling the behavior of cells on micro-patterned substrates have established the emer¬ 
gence of coherent rotation of cells under confined conditions nans]. Similar to these studies, our model also shows 
the emergence of a persistent mode of rotation for a group of cells (N = 140) when confined on a circular substrate 
(k c = k t = 10, £ = 1, vq = 1, fi = 1) (Video SI). While the theory of active elastic systems attributes the onset 
of rotational motion to energy transfer to the lowest modes [51/80!, a systematic analysis of this phenomenon in 
the context of epithelial sheets remains to be performed. Using our model, we demonstrate that rotation is indeed 
the preferred mode of motion for tissues confined in circular geometries—this mode of CAM is very different than 
that observed in bacterial suspensions [20] (also see SI Text S2). Figure [2]] c) illustrates the quantification of this 
rotational motion in terms of mean vorticity of the system (See materials and methods). After an initial transient 
mode, cells start to rotate steadily as evidenced by the constant value of the mean vorticity of the system. The 
onset of rotation depends on the parameter £, which reflects the tendency of the cell’s polarization to orient along its 
velocity (Fig. Efa)). The greater the value of £, higher is the tendency of polarization vector to reorganise and align 
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along the velocity vector, resulting in faster initiation of coherent rotation of cells (Fig.[2fb)). Figure[2fd) emphasizes 
this by plotting the scalar product of polarization vector and velocity vector (p.v) as a function of time. From the 
figure it is seen that, as the value of £ increases, coordination between p and v is builds up faster resulting in a faster 
approach to steady state of motion. We would also like to emphasize that, for larger values of £, the time scale for 
polarization evolution can be faster than the relaxation of a few long wavelength radial modes (see SI Section Text S2 
and SI Video S3). In this case, some long wavelength radial modes can be sustained during the coherent rotation 
and the tissue can exhibit radial movements that are similar to those observed by Deforet et al. [83]. Additionally, 
as the confinement radius R for the tissue increases, these radial movements become prominent even at lower values 
off ( Section Text S2 of SI and SI Video S19). This is because, larger the system size, lower is the stiffness of long 
wavelength radial modes, and hence slower is their decay. This behavior of increasing radial velocity for the tissue 
with increasing confinement size is also observed by Deforet et al. in their experiments (see SI Fig. 4 of Ref. [83]). 
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(a) Time evolution of p and v for £ = 0.1 
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Figure 2. Coherent rotation of cells on circular geometry, (a) The time evolution of polarization vector, p 
and velocity vector v is shown for £ = 0.1. The evolution rule for polarization is chosen in such a way that, from an 
initial random orientation, p will try to orient along velocity vector with time, (b) The coordination between p and 
v is decided by the parameter £. The higher the value of £, higher is the tendency of p to orient along v. The 
orientation of p and v at steady state for £ = 0.5 and £ = 1 are also shown, (c) Mean vorticity for systems with 
different £ is plotted as a function of time, (d) The tendency of polarization vector to orient with velocity vector is 
shown by the plot between p.v and time. As the value of £ increases, value of p.v approaches 1, indicating perfect 
alignment between two vectors, (e) A plot of velocity corFelation length for varying system size shows that 
correlation length equal to the confinement size, (f) A plot of correlation function with time shows that the velocity 
correlation length increases with time, till the coherent rotation sets in. 





















It was reported by Doxzen et al that, for tissues with confinement size greater than the velocity correlation 
length (« 200 /i m), there was no onset of CAM within the observation window of around 48 hours [9]. However, we 
find from our simulations that irrespective of tissue size (R), the tissue always reaches the steady state of coherent 
rotation (see Fig. [2]). In other words, we find that the steady state velocity correlation length is set by the size of the 
confined tissue. However, the time required to reach the steady state is higher for larger tissues (see Figs. Efe)-(f) 
and SI Fig. S5). This increase in the time required to reach the steady state may be attributed to the presence of 
a greater number of long wavelength modes for the larger system, as described above. The presence of these modes 
would interfere with the transfer of cellular motility to the rotational mode. We can reconcile our simulation results 
with the experimental observations by noting that, as the time required for setting the coherent motion is greater 
for larger tissues, the tissue is likely to be perturbed by certain unknown factors (e.g., cell proliferation) in that 
additional time. The resulting mechanical and polarization perturbations may, therefore, further delay the onset of 
coherence with respect to the experimental time window, or make CAM infeasible. We predict that in the absence 
of perturbations, even a large confined tissue can undergo CAM. These predictions differ from the observation of 
finite velocity correlation lengths of around 10 cell lengths in unconfined tissues (e.g. Refs. mm), wherein different 
boundary conditions (e.g., leader cells, high cable tension, etc) are likely to lead to qualitatively different behavior 
from that of confined tissues. 

Collectively, these results illustrate the effect of confinement in inducing coherent angular motion. Under in vivo 
conditions, such confinement may be provided by non-motile cells [35] possessing higher substrate frictions than 
motile cells (see SI Text S4 and SI Figs. S3, S4, Video S4-S11). Under these conditions, the efficiency of coherent 
motion is dictated by the ratio of substrate frictions between the two cell types. 

Cell crowding leads to fluidisation of tissue 

As the presence of a rotational mode of migration under confinement is well established by now, we focused our 
attention in understanding the characteristics of that motion in detail. Studies by Doxzen et. al. have shown that 
the movement of small circular tissues under confinement is similar to solid body rotations with angular velocity uj 
equal to |j^, where R is the radius of circle [9]. Further, the linear relationship between velocity and radial distance 
for rotating cell collectives obtained by multiple research groups support the argument of solid body rotations [9 J93j- 
However, what factors influence this solid-like tissue behavior has not been addressed. Here, we show that cell 
density is one such parameter dictating the nature of tissue behavior. As shown in Video SI, at lower cell densities, 
system behaves as an elastic solid with negligible neighbor changes and a linear velocity versus radial distance 
relationship (Fig. [3fa)). Increase in number of cells in the system while keeping the size R constant, i.e., increase in 
cell density, leads to an interesting phenomena. Increase in cell density alters the nature of the velocity versus radial 
distance relationship and induces a transition from solid-like behavior (N = 140) to that like a fluid (N = 170). 
Specifically, with increase in cell density, the linear velocity versus radial distance curve becomes more saturating. 
At the highest cell density (N = 170), the velocity plateaued to vq = 1 at the edges. One of the probable reasons for 
this change is the large shear that the system experiences at such densities, as evident from the relative sliding of 
cells past each other (Video S2). Quantification of the shear strain rate (e xy ) from the rate of deformation tensor as 

e xy = \ was performed to obtain additional insight into the magnitudes of shear experienced by the cells 

at various cell densities. A plot showing the variation of principal shear strain rate as a function of radial distance 
shows that with increase in cell number, the shear in the system also increases (Fig. [3^b)). Collectively, the above 
numerical results indicate that the number density of cells alters the behavior of system; i.e, at lower cell densities, 
system behaves like an elastic solid and at higher cell densities, system becomes more fluid-like. 




(C) 




Number of cells (N) 


Number of cells (N) 


Figure 3. Cell crowding leads to fluidisation of tissue, (a) The relationship between velocity and radial 
distance is examined for varying number density. Keeping the values of other parameters same as in previous 
simulations, the absolute velocity, |v| averaged over time, after the system reaches steady state, is plotted as a 
function of radial distance for varying number of cells N = {140,150,160,170}. As the number density of system 
increases, the velocity-radial distance curve become less linear, indicating the presence of shear in the system, (b) 
Variation of principal shear strain rate along the radial distance plotted as a function of number density. Increase 
in shear rate with number density illustrates the fluidisation of tissue induced by cell density, (c) Vorticity of 
system decreases with increase in cell density, (d) Without considering the effect of contact inhibition, mean 
velocity of the system increases with number density. 


While studying the effect of cell crowding on the nature of coherent rotation, we assumed that the motile cell 
speed or the fraction of motile cells is not modified by cell density. Consequently, we find that the mean speed of 
the cells in the tissue increases with cell density (Fig. |3jd)). This finding follows from our observation in Fig. [3{a) 
wherein upon increase in cell density, the tissue fluidises, as a result of which more and more layers of the tissue 
move with speeds comparable to Vo = 1. On the other hand, when the tissue behaves elastically (for N = 140), the 
tissue rotates as a rigid body with cell speed comparable to vq at the edges, but significantly lower speed of cells in 
the interior. However, while studying the effect of cell density on velocity profile of the over-confluent tissue, the 
condition of contact inhibition observed experimentally |52] has not been taken into account. To mimic the condition 
of contact inhibition for a denser system, and reconcile the experimental observations of decrease in mean velocity 
with increase in number density [S], we have considered the following cases: (i) due to crowding, the self-propelled 
speed of cells can be smaller on account of cells forming smaller lamellipodia |9| (see Fig.[4{a)), or (ii) due to crowding, 
a fraction of cells are possibly not motile (see SI Fig. S6). Both of these effects are feasible due to contact inhibition 
of motility in crowded tissues. For both cases, as expected, we observed reduction in mean cell speeds. Additionally, 
we can also see from Fig. |4{b) that the tissue shows fluidisation for value of vq as low as 0.3; only at really low 
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vo = 0.1 does it recover back its elastic behavior. Thus, for appropriate values of vo at large TV, we can observe lower 
mean cell speeds, concurrently with a fluid-like behavior for the overall tissue. 


Figure 4. 

(a) Mean 



(b) 
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Cell motility dictates the fluidized behavior of tissue. 

velocity for varying values of cell motility (uo). (b) Normalized velocity-radial distance plot for 

values of vq for N = 170. 


varying 


Exact steady state solution when the tissue is a linear, isotropic, homogeneous, elastic 
solid 

It appears from the above findings that N and vo are two parameters that tune the solid to fluid transition in 
coherently rotating tissues. In order to gain insight into this transition, we first analytically obtain mechanical 
steady state of the tissue by modeling it as a homogeneous, linear, elastic solid. This continuum description seems 
reasonable when cell-cell connectivity in the tissue is maintained during coherent rotation m • For such a tissue 
constrained within a circular patch and with zero tangential tractions, we can obtain one particular steady state 
solution, based on the following assumptions: 

1. The solution is radially symmetric, i.e., there is no 0 dependence 

2. The polarization p is aligned in the 0 direction 

3. All cells are motile 


Due to ensuing radial symmetry, the confinement is expected to induce isotropic compression only, i.e., a r = ctq = p, 
which should not interfere with the shear strain (stress) in the tissue due to motile forces. 

The easiest solution to visualize is a small elastic displacement u r and uq superimposed on a rigid body rotation 
with angular speed uj. If this solution is indeed possible then, p and cell velocities will be both aligned in the 
tangential direction. The continuum form of polarization evolution equation (see Eq. SI) would be, 


Dp 

Dt 


= £(P x v).e 2 p_L 


(5) 


for the current model, where D/Dt represents the co-rotational material derivative for the elastic sheet [2H1IH3] . We 
look at the steady state solution when p would not vary temporally. 

The equation of equilibrium, respectively, in the radial and the tangential direction for the current elastic sheet 
with above conditions will be [87) : 


Eh f d 2 u r 1 du r u r \ Eh d fdu r u r \ 

2(1 + v) \ dr 2 r dr r 2 ) 2(1 — v) dr \ dr + r J 

Eh fd 2 ue 1 due uq \ vo / uor\ 

2(1 + v) \ dr 2 ^ r dr r 2 ) + /a s \ vq J 
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In the above set of equations E , v and h are, respectively, the Young’s modulus, Poisson’s ratio, and thickness for 
the sheet—the connection between these values and the parameters used in the simulations is discussed in SI Text SI. 
The parameters Vo and p s are the self-propelled speed and effective motility per unit area of the tissue. Since we 
presume that all cells are motile, vq is the essentially same as the self-propelled motility value used in our simulations 
for the tissue. The parameter p s is related to the motility of single cell as n s p = //, where p is the cell density, or 
number of cells per unit area of the tissue. 

The angular velocity uj is an unknown in this problem and can be obtained as follows. The equation of equilibrium 
for the tissue in the simplest form is: 


V.<7 + ~(vo 


— ujr)t = 0. 


( 8 ) 


Taking a cross product on both sides with r, the position vector with respect to the center, and integrating this over 
the entire area of the circle we get: 


e z 



V.crdA + e z . 



ujr)r x tdA = 0. 


(9) 


Since the tangential traction is zero, by design, on the boundary, the first term of this equation reduces to zero by 
the divergence theorem. The second term can be simplified, further, due to the presumed radial symmetry to give 
the following expression for u ): 


uj = 


4vo 
3 R ’ 


( 10 ) 


where R is the radius of the confined tissue. This derivation is similar in essence to that done in Ref. [9]. Substituting 
this value for u) we can now solve the two equations subject to two boundary conditions: 


u r (R) = 0 (confinement), 


d_ fu e \ 
dr ' 


-) 


R 


= 0 (shear traction). 


(ii) 


Using these boundary conditions, and noting that, by symmetry iz r (0) = uq(0) = 0, we get the following solution for 
the displacements (with respect to the undeformed configuration) 


u r = 0 

V 0 p{l + v)R 2 f r \ 2 / r \ 

Uf> ~ 3Ehp \RJ \R J 

The internal shear r r e (per unit height h) for the tissue is given as 

Eh f due_ __ uq\ _ vqR_ /_r \ fr_ _ 
Tr0 2(1 + z/) \ dr r ) 3p s \RJ \R 

The maximum value r max of r r e happens at r = R/2 and given as 

_v 0 pR ^ _v 0 pR(l + v) 

tnax — ana e max — , 

Yip llpEh 


( 12 ) 


(13) 


(14) 


This simple expression, in combination with simulations, gives us significant insights into the variables and mecha¬ 
nisms that influence the behavior of the coherently rotating tissue. 

It can be seen from Eq. [IT] that the maximum shear strain in the tissue, assuming it to be linear, elastic, 
homogeneous, is directly proportional to both vq and p ~ N. Because shear strain governs transition from solid to 
plastic flow m , we expect the quantity Vo x N to govern the transition of the tissue from solid-like to fluid-like. 
Since we observed solid-like coherent rotation of the tissue for N = 140 and vo = 1, we can expect the tissue to 
behave in a similar manner for N = 170, if the cell motile speed is taken as vq = 140/170 ~ 0.82. It can, however, 
be seen from Fig. 0] (b), that even for values of as low as 0.3, the tissue still undergoes fluidisation. Only when 
reduces to values lower than 0.1, does the tissue recover back its solid-like behavior. This implies that there is 
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possibly a density dependent shear threshold which controls the solid to fluid transition of the tissue. This can be 
achieved if confinement introduces a density dependent shear pre-strain in the tissue, thus apparently altering the 
critical threshold. 

The underlying mechanism of tissue fluidisation at higher densities can perhaps be understood from analysing the 
deformation pattern of cell triangles in the tissue. In our study, the cells, represented by their centers, are confined 
in a circle of a given radius R. The presence of confinement results in pre-straining of cell-cell connections (springs). 
Unlike in a homogeneous material, due to discrete nature of this system, pre-straining by circular confinement is not 
uniform; instead, it leads to non-uniform deformation of the springs resulting in the presence of shear pre-strain in the 
system, which, by definition, implies distortion of the tissue. This distortion is further enhanced by the shear strain 
induced by the motile forces on the cells. Delaunay triangulation, which is used in our model to obtain/update the 
connectivity of cells, seeks to minimise distortion (shear) in the connecting triangles that form the tissue. Everything 
else remaining the same, crowding increases the amount of pre-strain, and hence the initial shear strain in the tissue. 
Thus a crowded tissue is more susceptible to connectivity update via neighbor changes (i.e., T1 transition), which is 
reflected as non-zero shear strain rate or fluidisation (see SI Sections Text S5 and Text S6 for detailed analysis). 

It may be noted that, in order to provide a general and more realistic, continuum description of the tissue here, 
one may need a more sophisticated model m • This is beyond the scope of this paper, due to the difficulty in 
both obtaining appropriate rheology that is compatible with the discrete model, and obtaining an analytical solution. 
We hence, present a simple semi-analytical case of a simple Newtonian fluid to demonstrate coherent rotation for 
a fluidised tissue, and resort to the simulation results to make any contact with experiments (see SI Text S3 for 
derivation). The analytical predictions of mean vorticity (u; = ^ in the case of solid-like and uj = ^ in the case of 
fluid-like (calculation shown in SI Text S3)) closely matched the simulation values, and predict a reduction in mean 
vorticity with increase in cell density (Fig. [3jc)). 

Effect of tissue size, cell stiffness and cell cohesivity on tissue fluidisation 

As seen from the previous section, in addition to providing rotational velocity, the continuum modeling also gives us 
a simple expression for maximum shear strain (stress) in the tissue (Eq. [14]). This equation gives us further insights 
into the possible behavior of the tissue. For example, this expression predicts that a tissue with larger R has greater 
shear strain, and is hence more susceptible to cross over the critical strain threshold and exhibit fluidisation. To test 
this prediction, we performed simulations with increasing i?, such that the number density of cells in the tissue was 
very close to the number density for the case R = 5,7V =» 130, where the tissue rotates as a solid. It can be seen 
from Fig. [5](a) and SI Video S19 that, though there is no fluidisation for R = 5, for larger i?, the tissue behaves in 
an increasing fluid-like manner—more and more layers of tissue were observed to move with velocity close to = 1. 
Thus, tissue can undergo fluidisation solely due to the influence of system size. The relatively larger values of cell 
speeds at lower radial distance is due to radial movement of cells (see SI Video S19), and is possibly related to the 
dominance of radial modes with increasing system size (SI Section Text S2). Thus, even though Eq. [14] does not 
exactly capture the tissue behavior with increasing system size, it provides us with pointers in the right direction, 
and concurrently exposes the shortcoming of describing the tissue as a solid-like material [9]. 
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Figure 5. Tissue size, cell stiffness and cell cohesivity influence the fluid-like behavior of tissue, (a) 

The relationship between velocity and radial distance is examined for three systems with varying radius, while 
keeping the number density approximately same for all. The number of cells in the systems are taken as N = {1170, 
520, 130} for R = {15, 10, 5}, respectively. The values of other parameters are chosen as that of the previous 
simulations. It is observed that, while keeping the number density constant, with increase in system size, the 
velocity versus radial-distance profile become less linear as more number of cells tend to move with a velocity 
comparable to this shows the presence of shear strain rate in the system ( see SI Video S19). (b) Increase in cell 
stiffness by increasing the value of compressive stiffness (k c ) of a system will make the system stiff and resulting 
rotational behavior will be more like a solid, (c) Reduction in cell cohesivity (k t ) leads to fluid-like tissue behavior. 

It can be noted from Eq. [14| that, the shear strain is, as expected, inversely proportional to the tissue stiffness. 
This implies that tissues with stiffer cells (k c ) and greater cell-cell cohesivity ( k t ) are less susceptible to cross over 
the critical strain threshold and more likely to exhibit solid-like behavior; the inverse would apply for tissues with 
softer cells. For the case R = 5,1V = 170, increasing the stiffness k c for a tissue from 10 to 100 results in a transition 
from fluid-like to solid-like coherent rotation of the tissue (Fig. 0[b)). Similarly decreasing the value of cell cohesivity 
(kt) also leads to fluid-like behavior of tissue. We can see from Fig. 0(c) that for N = 140 when k t = k c = 10, then 
the velocity profile being linear is an indication of rigid body rotation (as per the analytical solution for elastic solids 
shown in previous section). However, when k t is decreased from 10 to 1 while keeping k c = 10, then it is clearly seen 
that the tangential velocity as a function of radial position has saturating profile (as seen in the previous section 
for analytical solution for viscous fluid) indicating fluidisation. Thus the stiffness and cohesivity of tissue cells can 
independently control the nature of coherent rotation for the confined tissue. 

Coherent rotation in non-convex annular geometries 

While the above results of coherent rotation were obtained in circular geometries, it remains unclear if similar coherent 
rotation is also possible in non-convex geometries. Of the various non-convex geometries, annular rings are often 
observed in vivo in glands, ducts or tissues with lumen inside. Several studies have probed the collective behavior 
of cells inside annular geometries [SZIIMIES]- Since annular geometry represents the simplest non-convex geometry 
obtained from a circular shape, we next studied the coherence patterns in annular geometries and the influence of 
cell density. For this, simulations are done with outer and inner radius of annulus taken as 100 /im and 70 tim. 
Simulation with N = 100, k c — 10, k t = 10, £ = 1, v$ = 1, shows that here also, after a short initial transient mode, 
cells exhibit robust coherent rotation similar to that on circular geometries (Video S12). However, the pattern of 
coherent rotation is dictated by the stiffness of cell-cell adhesions. Specifically, for lower stiffness values (fc*), different 
cell layers in annular section may move in different directions (Video S13); in contrast, for higher stiffness values of 
cell-cell connections, cells move in a robust manner after an initial breathing mode (Video S12). 

Next, to test the effect of cell density on mean vorticity, simulations were performed on annular geometries with 
varying annular thickness, t and constant outer radius (R). For a constant number of cells (V), varying {- leads to 
change in number density, and is hence expected to influence the pattern of coherent motion. Consistent with this, 
distinct behavior is observed for two different values of N. Figure 0 shows the plot of mean vorticity of system as 
a function of {-for N = 100 (green curve) and N = 140 (blue curve); the other parameters are kept the same as 
in previous simulations. As seen from the plot, it is seen that with increase in the thickness of annulus, the mean 
vorticity of the system also increases, which matches with the findings of Li and Sun [93]. In addition to that, we 
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have also shown that for lower number of cells (N = 100), system behaves more like an elastic solid with minimum 
shearing between cell layers similar to that seen in the circular geometries. But it is interesting to note that for larger 
number of cells (N = 140), system behaves like an elastic solid at higher values. However, when the thickness of 
annular section decreases, cells become more compressed which leads to the fluidisation of system and as a result, for 
lower values, system behavior is more similar to viscous fluid. In order to have a better understanding, we have 
also calculated the analytical values of mean vorticity of system using the following equation: 

1 27r(i?out^out -^in^in) 

"=2 ,(^,-<) ' (15) 

from Stokes theorem m- Here R out and R[ n are the outer and inner radius and v out and are the outer and 
inner velocities, respectively. As seen from Fig.[6j the analytical values of vorticity (see earlier sections for analytical 
expressions for elastic and viscous calculation for v OVLt and ^ n ) closely follow the computed values and illustrate the 
dependence of vorticity on j- values. Taken together, our findings on circular as well as annular geometry imply that 
for different confinement geometries, cell behavior can vary between that of a perfectly elastic solid and a complex 
fluid depending on cell density. 



Analytical 

■ Elastic 

♦ Viscous 

Simulation 

—•— N = 100 

* N = 140 


Jl 

R 

Figure 6. Coherent rotation of cells confined in annular geometry. As observed for circular substrates, 
cells confined inside annular geometries also exhibit coherent rotation. Simulations done on an annular shaped 
geometry with outer radius, R and thickness, t show that mean vorticity of system decreases with increase in 
number density as in the case of a circle. Furthermore, simulations done with two sets of cell numbers (N=100, 
represented by green curve and N =140, represented by blue curve) show that at lower densities, system behaves 
like an elastic solid, roughly matching their values with analytical results (red points). For higher number of cells, 
cell behavior is more like an elastic solid for thicker sections. As the thickness of annulus reduced, cell state 
transitioned from an elastic solid to viscous fluid. The red and magenta curves showing the analytical values of 
elastic solid and viscous fluid respectively, are derived as explained in main text. 


Sensitivity of coherent rotation to cell division 

In addition to illustrating the role of confinement in inducing coherent motion, our results also demonstrate the 
critical influence of cell density in dictating the pattern of coherent motion. While cell density can be experimentally 
controlled in in vitro experiments, under in vivo conditions, cell density is controlled by cell division—a factor that 
was not taken into account in our simulations. While several computational studies have tried to understand how 
cell division influences morphogenesis EMU, the sensitivity of coherent motion to cell division remains unexplored. 
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Having demonstrated the robust influence of cell density in our simulations, we next probed the extent to which 
coherent motion is sensitive to changes in cell density effected by cell division events. Cell division can occur either 
synchronously (i.e., all cells divide at the same time) or asynchronously (i.e., cells divide at different times). During 
early stages of embryo development, cells generally exhibit multiple fast synchronized division, accompanied by a 
transition stage and subsequent slow non-synchronized divisions, with different cells having different stages of cell 
cycle HUBS]. In vitro the cell cycle of individual cells can be synchronized by serum starvation m • To study the 
sensitivity of coherent motion to cell division, we probed how changes in the total number of cells in a confined 
geometry would influence the pattern of rotation. For these studies, an annular geometry was chosen, as such 
geometries are biologically relevant |37j- Further, both synchronized and asynchronized cell division were introduced 
into an already rotating system to perturb the steady state rotational motion. 

A total number of 40 cells, which are below the level of confluence, were confined in an annular substrate of outer 
radius 100 /im and inner radius of 70 /im, and allowed to reach a state of coherent rotation (Fig.[7]^a)). Once this state 
was reached, cells were allowed to divide either synchronously or asynchronously. For implementing asynchronous 
cell division, each cell in the system was initially assigned a random cell cycle number between 0 and 1. In contrast, 
for synchronized division, the initial cell cycle number of all cells were set to 0. The time interval for cell division was 
assumed as 24 hrs for all our simulations, and represents the time taken by any cell to reach its cell cycle from 0 to 1. 
Any cell, on reaching a cell cycle number of 1, underwent division to form two new daughter cells, provided the cell 
area was above some critical area (see Materials and Methods). The two new daughter cells formed after division, 
were assigned equal and opposite polarization in random direction, and placed at ao/2 along the major principal axis 
of mother cell’s area (Fig. [Tfa)). For both synchronous and asynchronous division, cell division was stopped when 
the total cell number reached 80. 
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Figure 7. Synchronous cell division changes the sense of coherent rotation, (a) Cells are allowed to 
undergo division, where each mother cell on attaining a mature cell cycle number will divide and become two 
daughter cells with equal and opposite polarization (as indicated by black arrows inside the daughter cells). Two 
cases are analysed where cells are allowed to divide either synchronously or asynchronously. For both cases, initially 
starting with 40 number of cells, division is continued till number of cells become 80. (b) Even though incapable of 
change the direction of overall rotation, asynchronized cell division causes local perturbations in the pattern of 
rotation, which shortly dies down and system continues the steady rotational mode, indicated by the same sign of 
mean vorticity before and after cell division, (c) In the case of synchronous division of cells, the large perturbations 
introduced into the system in a short time is capable of inducing the change in the direction of rotation indicated 
by the opposite signs of mean vorticity before and after the division. Even though the reversal in the direction of 
rotation after synchronous division has not happened in all the cases analysed, we observed a preferential bias in 
the change in direction tendency. 


Interestingly, synchronous and asynchronous division were found to perturb coherent rotation to varying extents. 
Asynchronous division did not alter the direction of rotation, but only created some local disturbances, after which 
coherent rotation was fully established. This is clearly seen from the temporal profile of the mean vorticity ( Figs. 
[7|(a, b)) where transient fluctuations in the mean vorticity quickly die down and cells continue to rotate coherently. In 
contrast, for the case of synchronous cell division, on several occasions, the direction of coherent rotation underwent 
a change as observed from the change in mean vorticity values (Fig. [71(c)) . Statistical analysis revealed that out of 
200 independent simulations conducted, this reversal after synchronous cell division was observed for almost 60% of 
the cases, indicative of a preferential bias for the change in rotation (Video S14-S15). Together, these results suggest 
that while coherent rotation is insensitive to asynchronous division, synchronous division introduces a bias in the 
direction of rotation. However, the biological implication of this reversal remains to be established. 
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Effect of removal of confinement: cell stiffness and cell-cell cohesivity dictates invasion 
pattern from coherent motion 

Under in vivo conditions, the confinement assumed in our simulations, is generally provided by the surrounding 
extracellular matrix (ECM). For example, all epithelial tissues are surrounded by the basement membrane, which 
helps to maintain tissue organization and prevents cell invasion. However, the basement membrane is breached by 
epithelial cells which turn cancerous. Cancer cells are known to invade both as single cells and collectively mm- 
Since coherent rotation is sensitive to the properties of cell-cell contacts (i.e., k t and fc c , respectively) (Fig. [5j), we 
hypothesize that, the initial coherent rotation dictated by the properties of cell-cell adhesions has a distinct bearing 
on the eventual invasion pattern, when confinement is removed. To test this hypothesis, we have studied the invasion 
patterns formed when a coherently moving group of cells break their boundaries and invade to the surrounding matrix. 
For doing this, three conditions were chosen with the following combinations of k t and k c to mimic different properties 
of cells and cell-cell adhesions: k c = kt = 1 (i.e., soft), k c = 10, fct = 1 (i.e., medium stiff), and k c = 10, kt = 10 (i.e., 
stiff). The number of cells in each system was taken as 100 and the values of all other parameters were kept the 
same as that of other simulations. Once coherent rotation was set up in all the systems, the confinement was relaxed 
at t = 50 to allow for invasion. Consistent with our hypothesis, the combination of k c and k t were found to directly 
influence the nature of coherent motion (Fig. [8]^a-c) and Video S16-S18). For the soft and medium stiff systems, the 
extent of invasion (i.e., radial position as function of time) remained the same. However, contrary to the ‘soft’ case 
where cells scatter in all directions, for the ‘medium stiff case, cells move radially outward as clusters which remain 
connected. For the ‘stiff case, cells continue to rotate even after the removal of confinement. Together, these results 
demonstrate that the nature of coherent motion set by the extent of cell-cell cohesivity dictates the invasion pattern 
when confinement is removed. Also, the persistent rotation of stiff cells with stiff adhesions even after the removal of 
boundary shows that even though confinement is essential for the emergence of coherent rotation, depending upon 
the properties of the system, the presence of a confinement is not mandatory condition for the cells to continue in 
their coherent motion. 
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Figure 8. Cell stiffness and cohesivity dictate invasion pattern from coherent motion. Three different 
systems of cells are taken with different stiffness of cell-cell connections. Simulations for (a) a soft system with 
k c = 1 and k t = 1 (b) a medium stiff system with k c = 10 and k t = 1 (c) stiff system with k c — 10 and k t = 10. The 
number of cells in all the three cases are same and equal to 100. After reaching a steady state of rotation, 
confinement was removed at time, t — 50. The snapshots of cell migratory patterns at t = 55 and t = 60 are also 
shown. For the case of intermediate stiff system, cells migrate in clusters compared to softer system where cell 
invasion pattern is more scattered. At the highest stiffness, cells continue to rotate even after removal of boundary. 
The length scale for each set of figure is shown below them. 


Discussion 

Coherent rotation of cells is essential for a variety of physiological activities. In vitro studies have reported the 
occurrence of robust rotation of epithelial sheets under confinement. We developed a self propelled, cell center-based 
model to understand this phenomena. We have shown both numerically and analytically that, a two-way feedback 
between cell velocity and cell polarization, along with force transmission via cell-cell connections, is sufficient to 
produce persistent rotational mode of migration for cells under confined conditions. Though there are computational 
studies that demonstrate the possibility of CAM for tissues in confined geometries, we provide a logical and plausible 
explanation as to why this may happen. Similarly, we show the presence of additional hydrodynamic modes that 
deform the tissue in the radial direction, if the time-scale of relaxation of these modes is comparable to the time- 
scale of the orientation of polarization. We also predict that such radial modes are more likely to be observed for 
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larger size tissues since they can sustain low wavelength and slow decaying radial modes, whose stiffness is inversely 
proportional to the tissue size—this finding is consistent with the experimental observations of Deforet et al. [83] . 

Our model predicts that irrespective of the size of the confinement, the cells can reach the state of coherent 
rotation. This implies that the velocity correlation length for the tissue can be as large as the system size. However, 
Doxzen et al. reported absence of coherent rotation for larger confinement size and linked this observation with the 
system size being larger than the experimentally observed correlation length of ~ 10 cell lengths that was reported 
elsewhere [48]. This apparent contradiction regarding correlation length can possibly be explained as follows. Unlike 
confined tissues, the reported correlation length were obtained for systems having free boundaries. It is known that 
the presence of free boundaries leads to modifications in the tissue boundary conditions (e.g., leader cells, high cable 
tension, etc.) [68]. These modifications may influence the velocity correlation lengths for the tissue, and therefore 
lead to qualitatively different behaviors when compared with that of confined tissues. We would also like to point 
out that for the same type of cells, velocity correlation length can be influenced by mechanical perturbations—time- 
increasing velocity correlation length of up to 450 /im was reported in Ref. [69] for migration cells on deformable 
substrates. Thus the correlation length need not be an inherent property of a tissue type, and can be influenced 
by mechanical perturbations (such as confinement). In addition to these factors, perturbations in the form of cell 
division and cell death, can also influence the temporal dynamics of velocity correlation length of a confined tissue, 
and may be one of the reasons why coherent rotations were experimentally not observed for larger size tissues. 

Our results also demonstrate that a few experimentally controllable/observable variables—cell density (AT), motile 
speed (uo), system size (R), cell stiffness ( k c ) and cell cohesivity ( k t )—collectively tune the behavior of a coherently 
rotating tissue between that of an elastic solid and a complex fluid. We found that, in our model, tissue fluidisation is 
associated with increased shear strain in the system, which can be relieved by connectivity changes (T1 transitions). 
Such T1 transitions that are dependent on shear deformation of cells, have been also observed/modeled in a very 
recent work by Etournay et al. [96], thus providing biological relevance for our modeling. We found that upon increase 
in cell density the tissue showed an increasingly fluid-like behavior. This effect is not expected if the apparently 
radially symmetric confinement isotropically compresses the tissue. However, due to the discrete nature of the tissue, 
we observed that for greater cell density, the confinement induces larger distortion within the cells and makes the 
tissue more susceptible to fluidisation. Moreover, we also found that other factors such as increase in cell motility, 
increase in system size, and decrease in cell stiffness, which lead to increase in shear strains in the tissue, can make 
the tissue more vulnerable to fluidisation. Thus, our study identifies some of the potential parameters, which, as 
described above, could give rise to novel and hitherto unreported behavior for confined tissues. 

The basement membrane, which is found at the basal surface of epithelial cells is essential for tissue polarity, 
and maintains tissue structure by confining cells. In epithelial cancers, uncontrolled proliferation of cells leads to 
buildup of stress within the tissue. Subsequently, malignant cells breach the basement membrane and escape into the 
surrounding stroma. Cancer invasion through these matrices is dictated both by extrinsic factors (e.g., ECM density 
and organization) and by intrinsic factors. Among the intrinsic factors, our findings implicate cell division and 
functional nature of cell-cell contacts as two parameters influencing the coherent motion, and the invasion process. 
Our studies show that synchronous division introduces a bias in the direction of rotation, with a reversal in the 
direction of rotation observed in nearly 60% cases. Whether or not this reversal has any significance in invasion 
remains to be established. 

EMT, or epithelial to mesenchymal transition, refers to the complex process whereby immobile epithelial cells lose 
their cell-cell adhesions and get converted into motile mesenchymal cells HUES. EMT is relevant both to normal 
embryonic development and in carcinogenesis. During EMT, downregulation of the cell-cell adhesion protein E- 
cadherin is accompanied by upregulation of mesenchymal cadherins like N-cadherin which favour forming of transient 
contacts m- However, EMT is not an all-or-nothing phenomenon and cells can exist in partial EMT states. Such 
states have been reported in carcinosarcomas m ■ In contrast to EMT, cancer cells are also known to exhibit collective 
cell migration, where E-cadherin-positive cell-cell contacts are maintained m- It is likely that both cell scattering 
and collective cell invasion are outcomes of alterations in the physical behavior of cell-cell contacts. This was evident 
from the scattering patterns observed in our simulations when confinement was removed. When adhesions were soft, 
cells scattered in all directions with the formation and breakage of transient adhesions. This mode of invasion was 
closer to that of single cell invasion. However, when adhesions were medium stiff, cells scattered as uniform-sized 
clusters indicative of a collective mode of invasion. Our results thus suggest that the different modes of invasion 
observed in different contexts are dictated by the strength of cell-cell adhesions. When adhesions were very strong, 
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even upon removal of confinement, the cell layer expanded to release the confinement-induced compression, but 
continued to exhibit coherent rotation. Under this condition, activation of invasion as observed experimentally after 
removal of confinement 153 , is likely to involve mechano-chemical changes in the motility behavior of cells arising 
from the presence of free boundary. 

In conclusion, our framework of velocity-polarization coupling successfully recapitulates coherent motion in con¬ 
fined circular and annular geometries, demonstrates the influence of a few experimentally controllable variables— 
motile speed (i?o), cell density (TV), cell stiffness (k) and system size (R)~ in collectively dictating the pattern of 
coherent motion, and illustrates the effect of synchronous cell division on coherent motion. In addition, our model 
predicts the invasion patterns that arise due to coherent motion when confinement in removed. Future work can be 
focused on further improving the predictive power of our model by incorporating the effect of actomyosin contractility 
and substrate properties (e.g., stiffness). 


Materials and Methods 

Methodology 

For simulations, the scaling quantities for length, time and force are taken as ao, clq/vq and Vo/y respectively. Unless 
otherwise specified, values of ao, vq and fi are taken as 1 for all the simulations. The cells were represented by 
their centers, and their connectivity was obtained via Delaunay triangulation m, which produces least number of 
distorted triangles, i.e., triangles with least shear strain m • Delaunay triangulations are dual to Voronoi tessellations 
(see Fig. DJb)) and the Voronoi polygon for a given cell center can be modeled to represent the cell [89| l94j. It is 
unrealistic to obtain the areas of boundary cells directly from tessellation, since the Voronoi polygon for these cells 
can have a vertex at infinity m- To circumvent this problem, a row of dummy points are inserted at the boundaries, 
just to create well defined polygons for visualization, but they do not contribute to the dynamics of the system. Even 
though cells are connected to each other by cell-cell connections to form an apparently solid tissue, based on the 
dynamic position of the cells, this connectivity is constantly updated using Delaunay triangulation and may result 
in cell neighbor changes within the tissue. These modification of neighbors can be interpreted as the so called T1 
transitions, the specialised terminology for neighbor exchange in the context of foams and epithelia [951 (also see SI 
Section VI). Since confinement is experimentally shown to be essential for setting up coherent rotation, we model 
the soft confinement at boundaries by providing resistance of stiffness 3 k c at the edges, which will apply force on 
any cell trying to cross the boundary and thus prevent them from escaping. 

To begin with, cell centers were randomly distributed inside a confined zone of given dimensions, and were 
allowed to equilibrate, such that the velocities of all the cells was near zero—the cell-cell connectivity was obtained 
by Delaunay triangulation, as described. Since we model cells as self propelled particles, they were assigned a 
uniform motility, vo in random directions of their polarization after the equilibration stage. This led to the evolution 
of position and polarization of cells thereby setting up of the dynamics of the system. After a short initial transient 
state with random motility, cells started to rotate coherently and reached a steady state of motion. However, it is 
to be noted that the current formulation does not account for the effect of any noise in the system. For solving the 
set of differential equations numerically, we adopted forward Euler scheme that was implemented in Mat lab. After 
performing a detailed convergence study, a time-step of At = 0.001 was used for all the calculations. In order to 
quantify the angular motion of tissues, we calculated the mean vorticity of system derived from the antisymmetric 
part of velocity gradient matrix. Mean vorticity can be defined as where ^ = \ ~ §y)^ s the vorticity 

tensor, and u and v represent the velocity components in x and y directions at the time t. 

As already described, polarization of cells are initially randomly oriented. So it is logical to assume that there 
should not be any preferential bias in the direction of coherent rotation of cells. In order to verify this, a statistical 
analysis is carried out. Out of 100 independent simulations performed on both circular and annular geometry 
without cell division, almost equal number of clockwise and counter-clockwise rotations are obtained, which shows 
that there is no preferential bias in the system. Similarly, 100 independent simulations performed in an annular 
geometry with asynchronized cell division show no change in their direction of rotation after cell division. In order 
to check the statistical significance of the switch in the rotational direction on synchronous cell division, two sets of 
200 independent simulations are carried out. For the first set of simulations, the polarization of daughter cells are 
assigned equal and opposite in random direction, while for the second case, their polarizations are completely random 
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and independent. For both these cases, cells are observed to preferentially (around 60%) switch their direction of 
rotation after synchronous cell division, which indicate that the mechanical perturbations caused because of cell 
division may be the reason for the additional 10% bias in switch in direction. 


Supporting Information 

Text SI — Supplemental text showing the derivation of relation between spring constant 
(k) and Young’s modulus ( E ) for triangular network of springs. 

Let us consider a triangular network of springs with A, B, C as the initial position of cells and ao as the length of 
springs at equilibrium as shown in Fig. SI. Now, each cell is given a uniform stretch S so that their position has 
changed to A’, B’ and C’ respectively. For a network with six fold symmetry, the potential energy change associated 
with this stretching can be written as A F = 3 x \kd 2 [ 75 ] . where k is the spring constant. Assuming the system 
as elastic, isotropic and homogenous, the potential energy stored in the system as the result of deformation is equal 
to the strain energy of the system which can be written as V — \o%jtijV , where cqj is the stress tensor, e^- is the 
strain tensor and V is the volume of the element. Assuming the system as elastic and isotropic, the stress strain 
relationship can be written as eij = — where E denotes the Young’s modulus, v is the poison’s 

ratio and Sij is the unit tensor. For an isotropic stretch for a 2-D system, the normal strain components will be equal 
and shear will be zero. So e xx = e yy = e = ~ which will give the stress components as a xx = a yy = a = 

Now the total strain energy of the system will be U = (^~) 2 W Writing volume V in terms of the height of the 

element, h and equilibrium length of spring ao as V = y^(3)/4 x a^h, and equating strain energy with the change in 
potential energy, the spring constant will be k = 

Text S2 — Supplemental text showing the details of why coherent rotational motion 
is seen for a self-propelled elastic solid when the polarization vector for a cell has a 
tendency to align with the velocity of the cell. 

The equation of evolution for cell position and polarization, respectively, for the current system are (also see main 
text): 


v 0 Pi + pFi, 

€(Pi x Vi.e z )p± (SI) 

As per this evolution rule, the polarization of the cell has a tendency to align with its velocity. Additionally, the 
polarization direction also feeds into the velocity of the cell and tends to modify its speed and direction. Now, to 
understand, at least semi-analytically, the origin of the rotational motion under confinement for an elastic solid 
formed of self-propelling particles as given above, we use the following procedure motivated from the arguments 
in Refs. mm- We extend this reasoning to argue that even if the cells can exchange their neighbors, coherent 
rotational motion of the cell sheet is the most likely outcome. 

The elastic system of springs under circular confinement is free to undergo rigid body rotation about its center. 
Its translational degrees of freedom are, however, curtailed due to the confinement. Using ideas from structural 
mechanics, let us denote the stiffness matrix of the system in its rest (or stress-free) conformation as [K] |82|. The 
matrix [if], is symmetric, positive semi-definite, and has dimensions of 27V ce n s x 27V ce n s . If the displacement of the 
cells (nodes) from their rest position, in terms of column matrix of dimension 27V ce n s x 1, is {u}, then the equation 
of motion for the system in the complete matrix form can be written as: 

=v 0 {p}~ fJ.[K]{u}. (S2) 

Additionally, in this case, the polarization of each cell prefers to align with its velocity (see Eq. E3J. Since the stiffness 
matrix [K] is symmetric and positive semi-definite, it has 27V ce n s orthogonal eigenvectors {^} and corresponding 


drj 

dt 

dpi 

dt 
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non-negative eigenvalues AThe only zero eigenvalue is the one corresponding to the rigid body rotation mode 
{0o}; all other eigenvalues are positive. The displacement {xi} and velocity {tx} can then re-written in the form of 
eigen-modes as 


2-ATcells 1 

w = X! and ( S3 ) 

2=0 

2-/Vcells 1 

W = E ( S4 ) 

2=0 

where {p} is the polarization of all the cells, combined in the form of a column vector of size 2iV ce n s x 1. Expressing 
Eq.[S2] using the eigenmodes, we get the following set of equations in terms of eigenmode amplitudes 

^-=v 0 (4> j ){p}-nX j a j , (S5) 

where (00 is the eigenvector written in row format. This equation can be re-written as: 

=v 0 /3j - iiXjaj, (S 6 ) 

where /3j = (0j){p|. It seems safe to presume that for a certain time interval r ~ l/£, where £ is the response rate 
for polarization (see Eq. 1ST]) , the polarization remains almost constant. In this case Eq. [S5]can be solved to provide 
us the following solution for the amplitude oq of any mode j. 

otj = Y^-[l - exp(-fxXjt)], (S7) 

A jM 

and the corresponding velocity of the mode is given by 

dij — vo/3j exp(—/xAjt), (S 8 ) 

where we have assumed zero initial conditions for ctj. It is very clear from Eq. [S 8 ] that modes with larger A, i.e., 
greater stiffness or small wavelength, would decay faster as compared with the modes with smaller A. The smallest 
A possible for the current system is Ao = 0, corresponding to the mode that involves rigid body rotation of the 
tissue. This implies that the amplitude a o corresponding to pure rotation, and more importantly the angular speed 
do (vqPo) increases with time. The motility parameter p sets the rate at which the energy is dissipated from mode 
j 7 ^ 0, whereas the polarization orientation constant £ sets the rate at which the energy is pumped into each mode 
(in the form of /3j). 

Let us first look at the case with medium value of £ « 1, i.e., £ « (i\\. The polarization column vector {p} for 
the system can also be written in terms of eigen-modes as follows: 

2-^Vcells 1 

{p}= E M<hh (S9) 

2=0 

As can be seen from Eq. IS 8 l for the case of £ ~ /xAi, the velocity component corresponding to 0o increases, whereas 
the components corresponding to other modes essentially decay to zero—in the very least they do not grow as fast 
as ao- If may be noted that, since, the polarization vector for each cell has unit magnitude, the consolidated column 
vector will additionally satisfy 


2 -ZV"cells 1 

<P>{P}= E A 2 = iV cellB. (S10) 

2=0 

Hence, the fact that do is the dominant mode will ensure that 0o will increase to some bounded steady state value 
that would depend on the polarization orientation parameter £, since, as per our polarization rule the polarization 
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of a cell tends to align with its velocity. It may, however, be noted that since the polarization of each cell is a unit 
vector, in addition to a dominant /?o, some other components would also remain non-zero. This means that, as per 
Eq. [SOl some energy will keep getting pumped in a few other modes i. Nonetheless, do will be the only component 
that would increase steadily as per Eq. [S8]— other components would decay. 

We now examine two extreme limits of polarization orientation constant £. When £ is very small<< 1), the 
response of p to velocity v is slow. As a result, the polarizations of the cells would lag behind in their bid for 
orienting with the velocity (see Fig. 2(a), (b), (d) of the main paper), resulting in a smaller steady state value for 
/3q. However, as described above, even a small component /3 q of the polarization field would be sufficient to sustain 
steady do, and hence rotation—the other modes dq would not be sustained despite having non-zero f3j values in 
some of the modes. The angular velocity of the tissue (uj ~ voflo) would be, of course, small as is seen in Fig. 2(c). 
In the other extreme limit when £ 1, the rate at which energy is pumped in the modes is faster than the rate at 

which it is typically dissipated for the j th mode (^ fi\j ). In this case, we can see (Video S3), a perfect rotation of the 
tissue about the center is not obtained—the centre of rotation keeps changing constantly, confirming that {</>o} is not 
the only mode that is invoked. Indeed, as can be seen from Video S3, a few radial modes are also excited, and are 
reminiscent of such movements observed in Ref. [83]. Our calculation predicts that, if the cells are highly cohesive, 
in which case their polarization can evolve faster [M] (£ 1), we are expected to see these non-rotational, radial, 

modes. Similarly, since the stiffness of the long wavelength modes (A j) is inversely proportional to the system size 
(e.g., in 1 — D, Xj ~ j 2 /L ), we can see such modes for larger system even if £ « 1. In fact, such movements are also 
reported in the experiments of Deforet et al [83] and the simulations of Ref. [93] for larger system sizes—the authors 
attribute these modes to the lack of the strength of persistent force (vo//jL in our case). Our calculation, however, 
gives a clearer understanding for the origin of these movements. 

The previous argument hinged on the tissue having a well defined stiffness matrix [K], which in turn depends on 
having a system of cells with fixed connectivity. However, in our model we allow for the cells to change their neighbors 
and release internal stress. If that happens, the stiffness matrix of the tissue gets modified to [if'], depending on 
the rate at which the cells change their neighbors. As a result, the eigenvectors of the system now get modified to 
{</>'}. The only eigenvector that is, however, most certainly common to the two systems is the one corresponding 
to the rigid body rotation {0 O } = {^o}- Consequently, though there will be perturbations to in the form of new 
/3', the steady pumping of energy to the rotational mode will continue. The system is, hence, expected to achieve 
rotational coherence even if the cells (nodes) are allowed to change neighbors. Though this argument is not rigorous, 
it is consistent with the results of our simulations, and indeed seems plausible. 

There are a couple of things that we did not account in the above derivation: (i) pre-stress in the system due 
to crowding, and (ii) finite rotation effects. The pre-stress effect is too complex to be accounted for in this simple 
derivation—we leave it for future work. The effect of finite rotation seems to be a secondary effect. We see (Video SI) 
that when the cell sheet, by and large, behave elastically the coherent rotation is initiated before any finite rotation 
actually happens in the system. As a result, we think that the finite rotation effect is secondary, and if required can 
be incorporated by moving to a co-rotational frame of reference (as is done in section deriving the analytical solution 
for elastic solids in main text)—it should not affect essential mechanics of the problem. 

Text S3 — Supplemental text showing the details of exact steady state solution when 
the tissue is a viscous fluid. 

Since, in our model, the cells are allowed to change neighbours and relax the internal stress, then depending on the 
internal strain/stress, the tissue can indeed behave more like a fluid than like an elastic solid as described in the 
main paper. We hence, present a simple semi-analytical case of a Newtonian fluid to demonstrate coherent rotation 
for a fluidised tissue, and resort to the simulation results to make any contact with experiments. A more realistic, 
description of the tissue as a complex fluid is beyond the scope of this paper, due to the difficulty in both, using an 
appropriate rheological model, as well as in obtaining analytical solutions. 

In this case, we seek to obtain a radially symmetric solution such that both the polarisation p and velocities 
are aligned along the tangential direction (i.e., v r = 0). To do so, we write a very simple form for the equation of 
equilibrium as is given below. The equation for polarisation evolution remains the same as before (Eq. SI), and if 
we can find a such a solution, then we have found one possible steady state solution. 
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The constitutive equation for the epithelial sheet that is modelled as a viscous fluid is written as 

a = T] (Vv + (Vv) T - (V • v)l) + r] v (y • v)l, (Sll) 

where rj is the 2 — D shear viscosity, and r\ v = 3?? in 2—D [86] . The equation of equilibrium in radial direction will 
be trivially reduced to zero for the solution that we are looking for. In the 6 direction, the equation for tangential 
velocity vq will become, 


2 V 




V0_\ 
r 2 ) 


H-(vo - vo) = 0, 

f^s 


with the boundary conditions 


vq (0) = 0, (symmetry) 


f(-) 

or V r /r 


0 (shear traction at boundary). 


(S12) 


(S13) 


This equation seems to have a complicated closed-form solution in terms of Bessel and Hypergeometric functions. 
However, we can solve this problem numerically. To do that we will first non-dimensionalize the equation: all 
velocities are expressed in terms of vq and all lengths in terms of R. The equation then simplifies to 


1 d ( dv s \ Vg 
r dr V dr ) r 2 + 


vg) = 0, 



(S14) 


The quantity y /JRjj = Rh is the hydrodynamic length [86] , and the ratio a is the relative size of confinement disc (R) 
with respect to the hydrodynamic length R The non-dimensional Eq. lS14l equation can be easily solved numerically. 
For low values of a the solution seems to be very similar to a rigid body rotation. When a becomes larger, the velocity 
initially increases with r and then saturates to vq. One such plot for vq with respect to r is shown in Fig. S2. 


Text S4 — Supplemental text showing that passive confinement is also capable of in¬ 
ducing coherent rotation. 

In order to test whether cells can coherently rotate under passive confinements instead of geometrical confinements, 
simulations were performed with three initial positions of the ‘active’ cells: in the center of the circular pattern 
(Figs. S4 (a)-(c), Video S4-S5), in the periphery of the circular pattern (Figs. S4 (d)-(f), Video S6-S7) and in an 
annular geometry (Fig. S3, Video S8-S11). Consistent with our hypothesis, our simulations with varying mobility 
ratios of active and passive cells suggests that differences in physical properties of cells (i.e., motility and mobility) 
can indeed induce coherent angular motion (Fig. S3, Video S8 - Sll). When the frictional properties of active and 
inactive cells were comparable (i.e., //active///inactive) = 1, active cells were found to intercalate into the tissue and not 
exhibit any coherent motion. However, a 100-fold decrease in mobility of passive cell induced segregation of active 
cells due to inability of active cells to move the passive cells, and led to onset of coherent motion. This was even more 
clear for //active//inactive = 1000 where the majority of the active cells remained stuck in their initial positions, i.e., 
either centrally (Fig. S4(c)), or peripherally (Fig. S4(f)) or along the annulus (Fig. S3(d)), and exhibited coherent 
angular motion. Together, these results indicate that passive confinement effected by differences in cell mechanical 
properties is sufficient to induce coherent rotation. 


Text S5 - Supplemental text showing a detailed analysis to understand mechanisms 
that govern fluidisation of the tissue during coherent angular movement. 

We had addressed various parameters and mechanisms that can lead to fluidisation of the tissue in Section “Cell 
crowding leads to fluidisation of tissue.” Here, we present a detailed analysis of those mechanisms. Note that in the 
discussion below we interchangeably use distortion and shear strain; distortion implies change in shape, and shape 
deformation is the hallmark of shear strain m- 
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Confinement contributes to the distortion of the tissue 

In the absence of any compression, the approximate size of a single undeformed cell in a tissue, assuming the cell shape 
to be approximately circular, will be 7r(ao/2) 2 . As a result, the critical number of cells that can be approximately 
packed in a circle of radius R is N c « R 2 /(ao/2) 2 . Using the canonical values for our model, R = 5 and ao = 1, 
this number N c « 100. Hence, in a loose sense, a tissue under this confinement with more than 100 cells is expected 
to have deformation. If the tissue were elastic and homogeneous, one would expect the confinement leading simply 
to uniform and isotropic compression [87] . However, due to its discrete nature, the tissue is susceptible to having 
shear deformations in addition. This is because, the circular confinement does not permit all the cells to have a 
preferred co-ordination number z — 6 [EE], thus leading to distortion of the tissue, since the triangles formed by 
cell-cell connections (springs) cannot be all equilateral. The presence of motility forces further distorts the tissue by 
contributing additional shear strains. In the steady state of coherent rotation, we hence, expect the distortion of 
cell-connection triangles to depend on cell density. We will see below that, these intuitive observations are indeed 
consistent with the findings of our simulations. 

As described in the main paper, the connectivity of cells is obtained by Delaunay triangulation [89]. We can plot 
histograms of the qualities of cell-connection triangles spanning the tissue in Fig. S7. The quality of a triangle is 
defined as ED]: 


Q 


4.V3A 

h\ + h\ + h 2 ’ 


(S15) 


where, and hs are the side lengths of any given triangle, and A is its area. The quality factor Q = 1 for an 

equilateral triangle (no distortion), and reduces with increasing distortion. It can be clearly seen from Fig. S7 (a) 
and (b) that the number of relatively distorted triangles Q < 0.9 (chosen arbitrarily) is larger for the case TV = 170 
when compared with the less crowded N = 140. The steady state shear strain is, hence, observed to be comparatively 
larger for the more crowded tissue. 

Delaunay triangulation updates the connectivity of cells if the updated triangulation has fewer “skinny,” i.e, 
distorted triangles m- Hence, in a statistical sense, a more distorted tissue has greater susceptibility for changing 
its connectivity, thus resulting in its fluidization. 

Since Delaunay triangulations are dual to Voronoi tessellations [91 (Fig. 1(b) of the main paper), cell distortions 
can be equivalently but better explained by looking at the statistics of Voronoi (cell) edge length of the tissue in 
Fig. S8 (a) and (b) for N = 140 and N = 170. Cells connected to each other share an edge. When the cells 
update their connectivity via Delaunay triangulation, this edge grows in the topologically orthogonal direction while 
transitioning through a point, or a four-way vertex (see Fig. S9). Hence, in a statistical sense, a tissue with a 
substantial number of small edges is more susceptible to T1 transitions (or neighbour change). Since, it can be 
clearly seen from Fig. S8 that, in its steady rotating state, the number of very small Voronoi edges e < 0.2 (chosen 
arbitrarily) for N = 170 case, is quite large when compared with the case N = 140, we have another cue as to why 
the tissue is more prone to fluidisation for N = 170, when compared with N = 140. 


Increase in cell-cell connection (spring) stiffness reduces distortion and makes the tissue 
resistant to fluidisation 

Upon increasing the stiffness of springs, though the initial shear pre-strain is not expected to change much due to 
the kinematic confinement provided by the circular boundary, the additional distortion caused by the motile forces 
would clearly be decreased. The decrease in overall distortion of the tissue for N = 170 can be clearly seen from 
Figs. S7(c) and S8(c). Due to the relative lower distortion, as compared to a softer tissue, the tissue is more resistant 
to neighbour changes and hence fluidisation. 


Neighbour changes via Delaunay triangulation happen in local patches in a time- 
sequential manner and releases local shear strain in the tissue 

For any configuration of cell centers, Delaunay triangulation, in essence, provides the cells with connectivity that 
produces the least number of distorted cell triangles. As demonstrated in Fig. S9 (also see SI Video S22), for 
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N = 170, the neighbour exchange happens at local “hot spots,” in a frequent, and time-sequential manner. A local 
T1 transition event release the local shear strains in the tissue, as and when the strain builds up, and a neighbour 
change can reduce the same. As also described in the SI section below, T1 transition happens in our model only 
if re-triangulation results in a tissue with lesser distortion. If we prevent the T1 transition from happening, by 
“locking” the connectivity of cells, i.e., by not updating the connectivity of the cells, then distortion starts building 
up at local patches (SI Fig. S10(b) for N = 170) and steadily increases. Updating the tissue connectivity by Delaunay 
triangulation, which leads to fluidisation of the tissue, prevents build-up of shear strain and corresponding stress 
in the tissue. Thus updating cell-cell connectivity via Delaunay triangulation is a mechanically very relevant way 
to release shear strains in a tissue. Note that “locking” the connectivity does not lead to building up of stress for 
N = 140 (SI Fig. S10(a) and SI Video S21) because any modification in connectivity would only increase distortion 
in the tissue. 


Skewness of cell-cell connection length distribution is a key quantity that decides the 
susceptibility of the tissue to fluidisation 

The description given in the previous section clearly demonstrates the correlation between the shear distortion 
caused in the tissue by confinement with fluidisation upon crowding in the tissue. However, this does not explain 
why fluidisation is seen for N = 170 and not for N = 140, though confinement related distortion is present in both 
cases (see Fig. S7). Below we clearly describe that though the presence of shear strain is necessary for T1 transitions 
in the tissue, in itself it is not sufficient. For T1 transition to happen it is important that the local distortion of 
cells reduces upon T1 transition. Indeed, for N = 140, there are triangles with distortion comparable for the case 
N = 170 (see Fig. S7). What then is the possible reason for the stability of N = 140, with respect to T1 transitions, 
despite having distortion in triangles (albeit fewer numbers) that is comparable to N = 170? The resolution to 
this question, again, lies with the degree of confinement. For N = 140, since the density of cells is relatively lower, 
in addition to short edges we expect to find a chunk of long connections with length « ao = 1 at steady rotation 
state. The presence of such relatively long edges amongst short edges would lead to triangles of the type shown 
in Fig. SI 1(b). Though these triangles indeed have shear distortion, connectivity update will lead to triangles of 
even higher distortion. When the number of cells is increased, the proportion of longer edges decreases, resulting in 
triangles shown in Fig. SI 1(b) that are not as robust to triangulation update. 

To quantify, the proportion of long connections to short connections in the tissue we obtain the quantified the 
skewness in the distribution of edges for N = 140,150,160,170 given by Pearson’s second coefficient Sk 2 given 
by [22]: 


Sk 2 = 3 


(mean — median) 
standard deviation 


(S16) 


As shown in Fig. SI 1(c), Sk 2 increases with cell crowding, i.e., greater the skewness in the connection distribution, 
greater the amount of fluidisation of the tissue. This correlation is borne out consistently in other situations where 
shear can happen in the absence of crowding. 


1. Consider the case, when R = 10 and N = 560, this corresponds to cell density equivalent to R = 5 and N = 140. 
So, contrary to our exception, we indeed observe fluidisation of the tissue during steady rotation. This is very 
likely because the maximum shear in the tissue is proportional to R (see Eq. 14 and Fig. 5(a) of the main 
paper), as a result of which, the distortion due to motile forces can push the tissue over triangulation stability. 
The same behaviour is observed for R = 15 and N = 1080 (Fig. 5(a)). The index Sk 2 is around 0.4 for both 
cases, which is greater than the value of Sk 2 = 0.3 for R = 5 and N = 140 (no fluidisation). 

2. If we increase the radius of confinement to R = 5.8 such that the cell density for N = 170 is similar to the case 
R = 5 and N = 140, we indeed observe that there is no fluidisation in the tissue. The value of Sk 2 reduces 
from 0.54 (N = 170, R = 5) to 0.16. 

3. If the stiffness of the cell connections is increased to k c = k t = 100, the fluidisation of the tissue is prevented 
and the value of Sk 2 reduces drastically from 0.54 to 0.04. On the other hand, if k t is reduced to 1, then even 
for N = 140 and R = 5 case, tissue undergoes fluidisation, and the value of Sk 2 increases from 0.3 for no 
fluidisation to 0.4. 
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To summarise the findings, we find that the skewness ratio of connection lengths is intricately linked with the 
susceptibility of tissue to fluidisation. 

Text S6 - Supplemental text showing the details of Delaunay Triangulation and T1 
transitions in our model. 

In our model, the cell connectivity is obtained via Delaunay triangulation from cell centre positions at every time 
step |93j[9l]. The connectivity of the cells is only modified when there is local shear distortion of the tissue, which 
can be typically relieved by a systematic exchange of neighbour pairs (Figs. S9 and S10). These modification 
of neighbours that can be interpreted as the so called T1 transitions, the specialised terminology for neighbour 
exchange in the context of foams and epithelia [95]. T1 transitions that are dependent on shear deformation of cells, 
have been also observed/modelled in a very recent work by Etournay et al. [96], thus providing biological pointer for 
distortion driven neighbor changes in a tissue. It may be noted that, in the case of Vertex Models (VM), that have 
been extensively used for modelling epithelial morphology during development, the T1 transitions typically happen 
when the edge-length shared by two-cells reduces beyond an arbitrarily chosen critical threshold, upon which the 
connectivity is updated E51- The existence of a small edge in a cell, naturally implies the presence of shear strain, and 
is quite in line with distortion criteria used by Delaunay triangulation (and the corresponding Voronoi tessellation). 
However, Delaunay triangulation has one advantage of updating the cell connectivity naturally without resorting to 
potentially ad-hoc assumption used by VM mm- 



Figure SI. Triangular network of cells 
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Figure S2. Numerical solution for tangential velocity for a = 1 and 100. When a = 1, i.e., Rh ~ R , the 

tissue rotates almost rigidly. On the other hand, when a = 100, i.e., when R^ Rh, then the velocity increases with 
r and then saturates to value close to vq. 
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Figure S3. Passive confinement inducing coherent rotation, (a) A tissue comprising of heterogeneous 
populations of cells with active (red) and passive cells (green) undergo coherent motion depending upon the 
relative properties of two populations. Active cells, which are initially embedded in the form of an annular ring in 
passive tissue try to intercalate into the tissue, if their mobility ratios are comparable. With decrease in passive cell 
mobility (increase in friction), active cells experience difficulty in penetrating out from their initial position which 
results in a robust coherent motion of active cells, (a) - (d) show heterogeneous tissues with different mobility 
ratios ranging from 1 to 1000. In order to estimate the rotational motion of cells in the heterogenous tissue, we 
choose mean vorticity as the quantifying parameter as in earlier cases. For lower mobility ratios, active cells try to 
drag passive cells also along with them and try to build up a total rotation of the system. As the mobility ratio 
increases, their ability to drag passive cells decreases and for highest mobility ratio, active cells exhibit persistent 
rotation with steady value of mean vorticity. The spreading characteristic of cells from their equilibrium position is 
measured using mean radial distance and center movement. 
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t = 0 t = 300 t = 300 

Figure S4. Different positions of active cells in passive tissue, (a) - (c) show cells placed at the center and 
(d) - (f) show the cells placed at the periphery of passive tissue. Simulations with two different mobility ratios (10 
and 1000) show that system behavior is the same as that obtained for annular positioning of active cells. As in 
previous case, here also we see that for lower mobility ratio, active cells try to intercalate into the passive tissue and 
drag the neighboring passive cells along them. Similarly, for higher mobility ratio, active cells do not penetrate 
much outside their initial position, instead they exhibit coherent rotation. 



Figure S5. Correlation length normalized with system size is plotted for increasing time. As the size 
of the tissue increases, the time taken to reach coherence also increases. 
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Fraction of motile cells(N f ) 


(b) 



Figure S6. Velocity profile when only a fraction of cells in tissue are motile, (a) With decrease in 
fraction of motile cells, the mean velocity of system reduces, (b) Velocity profile for varying fraction of motile cell 
shows that as the number of motile cells increases, system behavior will change from solid-like to fluid-like. 



Triangle quality factor Triangle quality factor Triangle quality factor 

Figure S7. Quality factor for cell triangles of a coherently rotating confined tissue for different 
physical conditions. All other parameters remaining the same, it can be clearly seen that, (a-b) the proportion of 
relatively more distorted triangles is larger for the denser, fluidised, tissue (N = 170) as compared to a less dense 
tissue (N = 140). (c) When the stiffness of cell-cell connections (springs) is increased, the number of distorted cell 
triangles in the tissue with N = 170 is significantly reduced. As described in the main paper, the tissue now 
coherently rotates like a solid. 


(a) (b) (c) 



Figure S8. Voronoi (cell) edge length distribution for different tissue parameters. All other parameters 
remaining the same, it can be clearly seen that, (a-b) the proportion of relatively small edge lengths is larger for 
the denser, fluidised, tissue (N = 170), as compared to a less dense tissue (N = 140). (c) When the stiffness of 
cell-cell connections (springs) is increased, the number of small edges in the tissue with N = 170 is significantly 
reduced. As described in the main paper, the tissue now coherently rotates like a solid. 
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Figure S9. During coherent rotation of a tissue with higher cell density (N = 170), the neighbor changes happen 
in local patches with relatively higher distortion that can be reduced upon local connectivity update via Delaunay 
triangulation. The corresponding Voronoi cells at these places are also shown. It can be seen that such patches 
keep appearing locally at different places in a time-sequential manner (also see the related SI Video S20). 
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Triangle Quality Factor Triangle Quality Factor 

Figure S10. (a) “Locking” the node connectivity for less dense tissue (.N = 140) does not modify the mechanical 
state of the coherently rotating tissue (see SI Video S21). (b) On the other hand, “locking” the connectivity of 
denser tissue (N = 170) leads to build up of distortion (red triangles) in the tissue. It can be seen from SI Video S22 
that, upon “releasing” the connectivity lock, many cells undergo neighbor changes to relieve their distortion (shear). 



Triangle edge length Triangle edge length Number of cells 

Figure Sll. (a-b) It can be seen that the distribution of cell-cell connection (spring) lengths, which is a measure 
of distance between cell centers, becomes more skewed towards lower spring lengths, as the number density of cells 
increases, (c) The skewness as a function of cell numbers is quantified using Pearson’s second skewness coefficient 
Sk 2 - When the connection length distribution is less skewed, it implies that the number of large and small springs 
are comparable to each other. This may indeed lead to the presence of some distorted triangles edges as shown, but 
a connectivity update would lead to greater increase in distortion and hence not performed. The reverse is true 
when the distribution is more skewed. 


Video SI. Emergence of coherent angular rotation of cells confined in a circular geometry. 

The parameters for the simulations are N = 140, k c = k t = 10, £ = 1, vq = 1. 
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Video S2. Higher number density of cells leads to the fluidization of tissue, indicated by the huge 
shear appearing in the system. 

The parameters for the simulations are AT = 170, k c = k t = 10, £ = 1, Vo = 1. 

Video S3. For higher values of £, center of rotation keeps on changing. 

The parameters for the simulations are N = 140, k c = k t = 10, £ = 10, vq = 1. 

Video S4. Active cells initially placed in circular pattern in the center of a passive tissue. 

The parameters for the simulations are AT == 154, k c = 5, kt = 1, £ = 1, vo (active cells)=1, mobility ratio = 10. 

Video S5. Active cells initially placed in circular pattern in the center of a passive tissue. 

The parameters for the simulations are AT = 154, & c = 5, &* = 1, £ = 1, Vo (active cells)=l, mobility ratio = 1000. 

Video S6. Active cells placed on the periphery of a passive tissue in a circular pattern. 

The parameters for the simulations are AT = 154, fc c = 5, &* = 1, £ = 1, Vo (active cells)=l, mobility ratio = 10. 

Video S7. Active cells placed on the periphery of a passive tissue in a circular pattern. 

The parameters for the simulations are AT = 154, k c = 5, kt = 1, £ = 1, vo (active cells)=1, mobility ratio = 1000. 

Video S8. Passive confinement inducing coherent rotation in system with active cells placed inside a 
passive tissue in an annular pattern. 

The parameters for the simulations are AT = 154, k c = 5, kt = 1, £ = 1, vo (active cells)=1, mobility ratio = 1. 

Video S9. Passive confinement inducing coherent rotation in system with active cells placed inside a 
passive tissue in an annular pattern. 

The parameters for the simulations are AT = 154, fc c = 5, &t = l,£ = l, v $(active cells)=l, mobility ratio = 10. 

Video S10. Passive confinement inducing coherent rotation in system with active cells placed inside 
a passive tissue in an annular pattern. 

The parameters for the simulations are N = 154, k c — 5, k t = 1, £ = 1, vq (active cells)=l, mobility ratio = 100. 

Video Sll. Passive confinement inducing coherent rotation in system with active cells placed inside 
a passive tissue in an annular pattern. 

The parameters for the simulations are AT = 154, k c = 5, k t = 1, £ = 1, vq (active cells)=1, mobility ratio = 1000. 


Video S12. Emergence of coherent motion in annular geometry. 

The parameters for the simulations are AT = 100, & c = 10, = 10, £ = 1, vq = 1. 


Video S13. The pattern of coherent motion is determined by the stiffness of cell -cell connection. 

The parameters for the simulations are AT = 100, k c = 1, k t = 1, £ = 1, no = 1. 


Video S14. Asynchronous rotation does not switch the direction of rotation. 

The parameters for the simulations are AT = 40, & c = 5, &* = !,£ = !, vq = 1. 
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Video S15. Synchronous rotation switches the direction of rotation. 

The parameters for the simulations are N = 40, k c = 5, k t = 1, £ = 1, vq = 1. 


Video S16. Breakage of boundaries for softer system. 

The parameters for the simulations are TV = 100, k c = 1, kt = 1, £ = 1, no = 1. 


Video S17. Breakage of boundaries for medium stiff system. 

The parameters for the simulations are TV = 100, k c = 10, k t = 1, £ = 1, t?o = 1. 


Video S18. Breakage of boundaries for medium stiff system. 

The parameters for the simulations are TV = 100, & c = 10, &* = 1, £ = 1, vq = 1. 


Video S19. Coherent rotation of cells in larger systems. 

As the confinement radius R for the tissue increases, radial movements become prominent. The parameters for the 
simulations are R = 15, TV = 1260, k c = 10, k t = 10, £ = 1, vq = 1. 


Video S20. T1 transition happening in a denser system at local patches. 

The parameters for the simulations are N = 170, k c = 10, k t = 10, £ = 1, vq = 1. 


Video S21. “Locking” the node connectivity for less dense tissue does not modify the mechanical 
state of the coherently rotating tissue. 

The parameters for the simulations are TV = 140, k c = 10, k t = 10, £ = 1, vq = 1. 


Video S22. “Locking” the connectivity of denser tissue (N = 170) leads to build up of distortion 
(red triangles) in the tissue. Upon “releasing” the connectivity lock, many cells undergo neighbor 
changes to relieve their distortion (shear). 

The parameters for the simulations are TV = 170, k c = 10, k t = 10, £ = 1, vq = 1. 
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